%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reproduces Appendix Figure A.3
%
% Cloyne, Dimsdale and Hürtgen (2024)
% "Are Tax Cuts Contractionary at the Zero Lower Bound?
%  Evidence from a Century of Data"
%
% James Cloyne, Nicholas Dimsdale and Patrick Hürtgen, 2024
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all
close all
clc

%% set path and solve model using dynare and Occbin

set(0,'DefaultLineLineWidth',2)
addpath('c:/dynare/5.5/matlab')

dynare Occbin_NK_extended console

%% generate draws from prior distribution

draw2=2000;         % number of actual draws
draws=3000;         % extra number of draws in case some draws are not valid

draw_i=1;
go_flag=0;
go = true;

[beta1,beta2]=betaAB(0.5,0.2);
[beta3,beta4]=betaAB(0.6,0.15);
[gam1,gam2]  =gammaAB(2,0.5);
[beta5,beta6]=betaAB(0.8,0.1);

while draw_i < draws

    display(['Draw ', num2str(draw_i)]);

    omegap=random('Beta',beta1,beta2);
    omegaw=random('Beta',beta1,beta2);
    chip=random('Beta',beta1,beta2);
    chiw=random('Beta',beta1,beta2);

    phipi=random('Normal',1.5,0.2);
    while phipi<1.05
        phipi=random('Normal',1.5,0.2);
    end

    phiy=random('Normal',0.125,0.05);
    while phiy<0
        phiy=random('Normal',0.125,0.05);
    end

    rhot1=random('Beta',beta5,beta6);
    gamma  = random('Normal',1.5,0.25);
    theta  = random('Beta',beta1,beta2);
    psi    = random('Beta',beta3,beta4);
    adj    = random('Normal',6,1.5);
    xi     = random('gamma',gam1,gam2);

    actualdraws.omegap(draw_i)=omegap;
    actualdraws.omegaw(draw_i)=omegaw;
    actualdraws.phipi(draw_i)=phipi;
    actualdraws.phiy(draw_i)=phiy;
    actualdraws.gamma(draw_i)=gamma;
    actualdraws.xi(draw_i)=xi;
    actualdraws.chip(draw_i)=chip;
    actualdraws.chiw(draw_i)=chiw;
    actualdraws.theta(draw_i)=theta;
    actualdraws.adj(draw_i)=adj;
    actualdraws.psi(draw_i)=psi;
    actualdraws.rhot1(draw_i)=rhot1;

    draw_i=draw_i+1;
end

tic

%% start priod predictive loop

draw_i=1;
draws_skipped=[];
iter2=1;
iter3=1;

while iter3 < draw2

    drop_draw=0;

    omegaw= actualdraws.omegaw(draw_i);
    chiw  = actualdraws.chiw(draw_i);
    omegap= actualdraws.omegap(draw_i);
    chip  = actualdraws.chip(draw_i);
    phiy  = actualdraws.phiy(draw_i);
    phipi = actualdraws.phipi(draw_i);
    xi    = actualdraws.xi(draw_i);
    gamma = actualdraws.gamma(draw_i);
    theta = actualdraws.theta(draw_i);
    psi  = actualdraws.psi(draw_i);
    adj  = actualdraws.adj(draw_i);
    rhot1= actualdraws.rhot1(draw_i);

    set_param_value('rhot1' ,rhot1);
    set_param_value('omegap' ,omegap);
    set_param_value('omegaw' ,omegaw);
    set_param_value('chip' ,chip);
    set_param_value('chiw' ,chiw);
    set_param_value('phipi' ,phipi);
    set_param_value('phiy' ,phiy);
    set_param_value('theta' ,theta);
    set_param_value('psi' ,psi);
    set_param_value('adj' ,adj);
    set_param_value('xi' ,xi);
    set_param_value('gamma' ,gamma);


    options_.occbin.simul.SHOCKS=shock_mat.matrix_1;
    opts_simul_.maxit=100;
    [oo_, out]  = occbin.solver(M_, oo_, options_);
    if out.error_flag %solution was not successful
        drop_draw=1;
    end

    options_.occbin.simul.SHOCKS=shock_mat.matrix_2;

    [oo3_, out]  = occbin.solver(M_, oo_, options_);
    if out.error_flag %solution was not successful
        drop_draw=1;
    end

    if drop_draw==0

        IRF.linear_Y=oo_.occbin.simul.linear(1:50+0,strmatch('y',M_.endo_names,'exact'))-oo3_.occbin.simul.linear(1:50+0,strmatch('y',M_.endo_names,'exact'));
        IRF.piecewise_Y=oo_.occbin.simul.piecewise(1:50+0,strmatch('y',M_.endo_names,'exact'))-oo3_.occbin.simul.piecewise(1:50+0,strmatch('y',M_.endo_names,'exact'));

        IRF.linear_R=oo_.occbin.simul.linear(1:50+0,strmatch('r',M_.endo_names,'exact'))-oo3_.occbin.simul.linear(1:50+0,strmatch('r',M_.endo_names,'exact'));
        IRF.piecewise_R=oo_.occbin.simul.piecewise(1:50+0,strmatch('r',M_.endo_names,'exact'))-oo3_.occbin.simul.piecewise(1:50+0,strmatch('r',M_.endo_names,'exact'));

        IRF.linear_RR=oo_.occbin.simul.linear(1:50+0,strmatch('rr',M_.endo_names,'exact'))-oo3_.occbin.simul.linear(1:50+0,strmatch('rr',M_.endo_names,'exact'));
        IRF.piecewise_RR=oo_.occbin.simul.piecewise(1:50+0,strmatch('rr',M_.endo_names,'exact'))-oo3_.occbin.simul.piecewise(1:50+0,strmatch('rr',M_.endo_names,'exact'));

        IRF.linear_TT=oo_.occbin.simul.linear(1:50+0,strmatch('tt',M_.endo_names,'exact'))-oo3_.occbin.simul.linear(1:50+0,strmatch('tt',M_.endo_names,'exact'));
        IRF.piecewise_TT=oo_.occbin.simul.piecewise(1:50+0,strmatch('tt',M_.endo_names,'exact'))-oo3_.occbin.simul.piecewise(1:50+0,strmatch('tt',M_.endo_names,'exact'));

        IRF.linear_INFC=oo_.occbin.simul.linear(1:50+0,strmatch('infc',M_.endo_names,'exact'))-oo3_.occbin.simul.linear(1:50+0,strmatch('infc',M_.endo_names,'exact'));
        IRF.piecewise_INFC=oo_.occbin.simul.piecewise(1:50+0,strmatch('infc',M_.endo_names,'exact'))-oo3_.occbin.simul.piecewise(1:50+0,strmatch('infc',M_.endo_names,'exact'));


        Simulation.y1(:,iter3) = 100*IRF.linear_Y;
        Simulation.y2(:,iter3) = 100*IRF.piecewise_Y;

        Simulation.r1(:,iter3) = 400*IRF.linear_R;
        Simulation.r2(:,iter3) = 400*IRF.piecewise_R;

        Simulation.rr1(:,iter3) = 400*IRF.linear_RR;
        Simulation.rr2(:,iter3) = 400*IRF.piecewise_RR;

        Simulation.tt1(:,iter3) = 100*IRF.linear_TT;
        Simulation.tt2(:,iter3) = 100*IRF.piecewise_TT;

        Simulation.invfc1(:,draw_i) = 100*IRF.linear_INFC;
        Simulation.invfc2(:,draw_i) = 100*IRF.piecewise_INFC;


        iter3=iter3+1;

    elseif drop_draw==1
        draws_skipped(iter2)=draw_i;
        iter2=iter2+1;
        drop_draw=0;
    end
    draw_i=draw_i + 1
end

toc

%% generate figures

xmin=0;
xmax=13;
stepsize=1;

start1=5;
end1=5+12;

h=gcf;
set(h,'Position',[50 50 1200 800]);
myplot=figure(1);
fanChart(0:size(Simulation.y1(start1:end1,:),1)-1, Simulation.y1(start1:end1,:)); hold on;
set(gca,'TickDir','in');
xlabel('Quarters')
plot([min(xlim()),max(xlim())],[0,0], 'k','LineWidth',0.25);
box off; grid off
set(gca,'XLIM', [-inf inf], 'gridlinestyle', '-.', 'gridalpha', 0.5, 'FontSize',18)
set(gca,'XLIM', [xmin xmax-1]);
set(gca,'YLIM', [-0.5 0.5]);
set(gca,'YTick',[-0.5:0.5: 0.5]);
set(gca,'XTick',[xmin:stepsize:xmax-1]);
myplot.PaperPositionMode = 'manual';
set(myplot,'Position',[50 50 1200 800]);
orient(myplot,'landscape')
filepath1 = fullfile('FigureA3_GDP_nonZLB');
print(filepath1, '-depsc'); 

h=gcf;
set(h,'Position',[50 50 1200 800]);
myplot=figure(2);
fanChart(0:size(Simulation.y2(start1:end1,:),1)-1, Simulation.y2(start1:end1,:)); hold on;
%title('GDP at ZLB')
set(gca,'TickDir','in');
xlabel('Quarters')
plot([min(xlim()),max(xlim())],[0,0], 'k','LineWidth',0.25);
box off; grid off
set(gca,'XLIM', [-inf inf], 'gridlinestyle', '-.', 'gridalpha', 0.5, 'FontSize',18)
set(gca,'XLIM', [xmin xmax-1]);
set(gca,'XTick',[xmin:stepsize:xmax-1]);
myplot.PaperPositionMode = 'manual';
set(myplot,'Position',[50 50 1200 800]);
orient(myplot,'landscape')
filepath1 = fullfile('FigureA3_GDP_ZLB');
print(filepath1, '-dpdf', '-fillpage');
print(filepath1, '-depsc'); 

h=gcf;
set(h,'Position',[50 50 1200 800]);
myplot=figure(3);
fanChart(0:size(Simulation.y1(start1:end1,:),1)-1, Simulation.r1(start1:end1,:));hold on;
%title('Nominal Rate away ZLB')
set(gca,'TickDir','in');
xlabel('Quarters')
plot([min(xlim()),max(xlim())],[0,0], 'k','LineWidth',0.25);
box off; grid off
set(gca,'XLIM', [-inf inf], 'gridlinestyle', '-.', 'gridalpha', 0.5, 'FontSize',18)
set(gca,'XLIM', [xmin xmax-1]);
set(gca,'YLIM', [-3 0.5]);
set(gca,'XTick',[xmin:stepsize:xmax-1]);
myplot.PaperPositionMode = 'manual';
set(myplot,'Position',[50 50 1200 800]);
orient(myplot,'landscape')
filepath1 = fullfile('FigureA3_BankRate_nonZLB');
print(filepath1, '-dpdf', '-fillpage');
print(filepath1, '-depsc'); 

h=gcf;
set(h,'Position',[50 50 1200 800]);
myplot=figure(4);
fanChart(0:size(Simulation.y1(start1:end1,:),1)-1, Simulation.r2(start1:end1,:));hold on;
%title('Nominal Rate at ZLB')
set(gca,'TickDir','in');
xlabel('Quarters')
plot([min(xlim()),max(xlim())],[0,0], 'k','LineWidth',0.25);
box off; grid off
set(gca,'XLIM', [-inf inf], 'gridlinestyle', '-.', 'gridalpha', 0.5, 'FontSize',18)
set(gca,'XLIM', [xmin xmax-1]);
set(gca,'XTick',[xmin:stepsize:xmax-1]);
set(gca,'YLIM', [-0.25 0.25]);
myplot.PaperPositionMode = 'manual';
set(myplot,'Position',[50 50 1200 800]);
orient(myplot,'landscape')
filepath1 = fullfile('FigureA3_BankRate_ZLB');
print(filepath1, '-dpdf', '-fillpage');
print(filepath1, '-depsc');

h=gcf;
set(h,'Position',[50 50 1200 800]);
myplot=figure(7);
fanChart(0:size(Simulation.invfc1(start1:end1,:),1)-1, Simulation.invfc1(start1:end1,:));hold on;
%title('Inflation at ZLB')
set(gca,'TickDir','in');
xlabel('Quarters')
plot([min(xlim()),max(xlim())],[0,0], 'k','LineWidth',0.25);
box off; grid off
set(gca,'XLIM', [-inf inf], 'gridlinestyle', '-.', 'gridalpha', 0.5, 'FontSize',18)
set(gca,'XLIM', [xmin xmax-1]);
set(gca,'XTick',[xmin:stepsize:xmax-1]);
myplot.PaperPositionMode = 'manual';
set(myplot,'Position',[50 50 1200 800]);
orient(myplot,'landscape')
filepath1 = fullfile('FigureA3_INFL_nonZLB');
print(filepath1, '-dpdf', '-fillpage');
print(filepath1, '-depsc');

h=gcf;
set(h,'Position',[50 50 1200 800]);
myplot=figure(8);
fanChart(0:size(Simulation.invfc1(start1:end1,:),1)-1, Simulation.invfc2(start1:end1,:));hold on;
%title('Inflation away ZLB')
set(gca,'TickDir','in');
xlabel('Quarters')
plot([min(xlim()),max(xlim())],[0,0], 'k','LineWidth',0.25);
box off; grid off
set(gca,'XLIM', [-inf inf], 'gridlinestyle', '-.', 'gridalpha', 0.5, 'FontSize',18)
set(gca,'XLIM', [xmin xmax-1]);
set(gca,'XTick',[xmin:stepsize:xmax-1]);
myplot.PaperPositionMode = 'manual';
set(myplot,'Position',[50 50 1200 800]);
orient(myplot,'landscape')
filepath1 = fullfile('FigureA3_INFL_ZLB');
print(filepath1, '-dpdf', '-fillpage');
print(filepath1, '-depsc'); 

h=gcf;
set(h,'Position',[50 50 1200 800]);
myplot=figure(11);
fanChart(0:size(Simulation.invfc1(start1:end1,:),1)-1, Simulation.rr1(start1:end1,:));hold on;
%title('Real rate at ZLB')
set(gca,'TickDir','in');
xlabel('Quarters')
plot([min(xlim()),max(xlim())],[0,0], 'k','LineWidth',0.25);
box off; grid off
set(gca,'XLIM', [-inf inf], 'gridlinestyle', '-.', 'gridalpha', 0.5, 'FontSize',18)
set(gca,'XLIM', [xmin xmax-1]);
set(gca,'XTick',[xmin:stepsize:xmax-1]);
set(gca,'YLIM', [-1 1]);
myplot.PaperPositionMode = 'manual';
set(myplot,'Position',[50 50 1200 800]);
orient(myplot,'landscape')
filepath1 = fullfile('FigureA3_RR_nonZLB');
print(filepath1, '-dpdf', '-fillpage');print(filepath1, '-depsc') 

h=gcf;
set(h,'Position',[50 50 1200 800]);
myplot=figure(12);
fanChart(0:size(Simulation.invfc1(start1:end1,:),1)-1, Simulation.rr2(start1:end1,:));hold on;
%title('Marginal tax rate away ZLB')
plot([min(xlim()),max(xlim())],[0,0], 'k','LineWidth',0.25);
set(gca,'TickDir','in');
xlabel('Quarters')
set(gca,'YLIM', [-5 25]);
box off; grid off
set(gca,'XLIM', [-inf inf], 'gridlinestyle', '-.', 'gridalpha', 0.5, 'FontSize',18)
set(gca,'XLIM', [xmin xmax-1]);
set(gca,'XTick',[xmin:stepsize:xmax-1]);
myplot.PaperPositionMode = 'manual';
set(myplot,'Position',[50 50 1200 800]);
orient(myplot,'landscape')
filepath1 = fullfile('FigureA3_RR_ZLB');
print(filepath1, '-dpdf', '-fillpage');
print(filepath1, '-depsc'); 

h=gcf;
set(h,'Position',[50 50 1200 800]);
myplot=figure(9);
fanChart(0:size(Simulation.invfc1(start1:end1,:),1)-1, Simulation.tt1(start1:end1,:));hold on;
%title('Marginal tax rate at ZLB')
set(gca,'TickDir','in');
xlabel('Quarters')
plot([min(xlim()),max(xlim())],[0,0], 'k','LineWidth',0.25);
box off; grid off
set(gca,'XLIM', [-inf inf], 'gridlinestyle', '-.', 'gridalpha', 0.5, 'FontSize',18)
set(gca,'XLIM', [xmin xmax-1]);
set(gca,'XTick',[xmin:stepsize:xmax-1]);
set(gca,'YLIM', [-1.5 0.5]);
myplot.PaperPositionMode = 'manual';
set(myplot,'Position',[50 50 1200 800]);
orient(myplot,'landscape')
filepath1 = fullfile('FigureA3_TL_nonZLB');
print(filepath1, '-dpdf', '-fillpage');
print(filepath1, '-depsc'); 

h=gcf;
set(h,'Position',[50 50 1200 800]);
myplot=figure(10);
fanChart(0:size(Simulation.invfc1(start1:end1,:),1)-1, Simulation.tt2(start1:end1,:));hold on;
%title('Marginal tax rate away ZLB')
plot([min(xlim()),max(xlim())],[0,0], 'k','LineWidth',0.25);
set(gca,'TickDir','in');
xlabel('Quarters')
set(gca,'YLIM', [-1.5 0.5]);
box off; grid off
set(gca,'XLIM', [-inf inf], 'gridlinestyle', '-.', 'gridalpha', 0.5, 'FontSize',18)
set(gca,'XLIM', [xmin xmax-1]);
set(gca,'XTick',[xmin:stepsize:xmax-1]);
myplot.PaperPositionMode = 'manual';
set(myplot,'Position',[50 50 1200 800]);
orient(myplot,'landscape')
filepath1 = fullfile('FigureA3_TL_ZLB');
print(filepath1, '-dpdf', '-fillpage');
print(filepath1, '-depsc');
